Contents

function hw4()

clear();
close all;
load('velb4degl.mat');

% Figure out how to partition this to get 50 DOF
dof = 50;
len_complete = size(velb4, 2);
points_per = floor(len_complete / dof);

for i=1:dof
    records.velocity(:,:,i) = velb4(:,(i-1)*points_per+1:i*points_per);
    records.time = 0:0.625:0.625*points_per; % times
    records.space = 0:2:2*219; % spacing
end


% Task 1: Present power spectral estimates of the horizontal velocity
% Ehat(sigma) at 50 degrees
% of freedom for two ranges separated by ?r=24m. Plot log-log on the same plot. Identify
% the swell peak, wind-wave peak. State the frequency in cycles/s, please.

% Get two records, 24m apart.
% From previous problems, we know that data in the spatial range 10:201 is
% good.  Each row is separated by 2m.  So, pick two rows 12 indexes apart.
% Index 10 = 20m (or maybe more... don't know what index 1 is, actually.)
% index 22 = 44m
%velocity_close(:,:) = records.velocity(10, :, :);
%velocity_far(:,:) = records.velocity(22, :, :);


% Get two records, 24m apart.
% From previous problems, we know that data in the spatial range 10:201 is
% good.  Each row is separated by 2m.  So, pick two rows 12 indexes apart.
% Index 100 = 200m (or maybe more... don't know what index 1 is, actually.)
% index 112 = 124m
% Both records will use the same frequency spacing, so we don't need
% different variables for that.
sample_interval = 0.625; % seconds
velocity_close = reshape(velb4(100,:), [], 1);
velocity_far = reshape(velb4(112,:), [], 1);
[ehat_close_real, ehat_close, freq_bins_real, freq_bins, delta_f] = ...
    varspec_est_multidof(velocity_close, sample_interval, dof);
[ehat_far_real, ehat_far, freq_bins_real, freq_bins, delta_f] = ...
    varspec_est_multidof(velocity_far, sample_interval, dof);

figure();
clf();
loglog(freq_bins_real, ehat_close_real);
hold on;
loglog(freq_bins_real, ehat_far_real, 'r--');
% I don't think we have a single axis label here.
%ylabel('Variance of velocity (m^2/s * hours)')
xlabel('Frequency (cycles/sec)');
legend('Close Range (200m)','Far Range (224m)');
title('Spectral density estimates using two ranges');
ylabel('Variance spectral density of velocity (cm/s) ^-^2');


% Task 2: Present co and quadrature spectral estimates at 50 d.o.f. Plot log (frequency,
% cps) vs. linear ordinate.

[cross_spec, freq_bins] = crossspectrum(velocity_close, velocity_far, sample_interval, dof);

figure()
% I
co = real(cross_spec(1:60));
semilogx(freq_bins, co);
hold on;
% Q
quad = imag(cross_spec(1:60))
semilogx(freq_bins, quad, 'r--');


% Task 3: Present estimates of coherence and phase (in degrees) at 50 d.o.f.
% On the phase plot, re-plot the phase estimate in each frequency band as an asterisk (*),
% but only in those frequency bands where the coherence is significant at the 67%
% confidence level.
coherence = abs(cross_spec(1:60)) ./ sqrt(ehat_close_real .* ehat_far_real);

% angle(z) is the same as atan2(imag(z), real(z))
phase_degrees = angle(cross_spec(1:60)) * (180/pi);

figure()
semilogx(freq_bins, abs(cross_spec(1:60)));
figure()
semilogx(freq_bins, phase_degrees);




end

function [cross_spec, freq_bins] = crossspectrum(data_a, data_b, sample_interval, dof)
    % We get 2 DOF from each set of data (thanks to real/imag components.
    num_sets = ceil(dof / 2);

    len_complete = size(data_a, 1);
    points_per = floor(len_complete / num_sets);

    delta_f = (1/sample_interval) / points_per;  % same as fs / N
    freq_bins = 0:delta_f:(points_per/2 * delta_f);

    % Get the variance spectral densitity for each subset.
    for i=1:num_sets
        start_idx = (i-1) * points_per + 1;
        end_idx = start_idx + points_per - 1;

        % normalize FFT by n^2
        a_a(:,i) = fft(data_a(start_idx:end_idx)) ./ points_per^2;
        a_b(:,i) = fft(data_b(start_idx:end_idx)) ./ points_per^2;

        crosses(i,:) = bsxfun(@times, conj(a_a(:, i)), a_b(:,i));
    end


    % Get the average values (average of columns)
    cross_spec = mean(crosses) ./ delta_f;

    % get rid of the first (mean) value.
    cross_spec = cross_spec(2:end);
    freq_bins = freq_bins(2:end);

end

Variance Spectral Density Functions

Take the given data and provide a variance density estimate The freq_bins array will be in reciprocal units to the units of sample_interval. ehat has units of (data_units^2) per (1/sample_interval_units)

function [ehat_real, ehat, freq_bins_real, freq_bins, delta_f] = ...
    varspec_est_real(data, sample_interval)

    % First, do the FFT to get an array of Fourier coefficients
    a = fft(data);

    % Figure out our frequency bins.
    N = length(data);
    delta_f = (1/sample_interval) / N;  % same as fs / N
    freq_bins_real = 0:delta_f:(N/2 * delta_f);
    freq_bins = 0:delta_f:(N * delta_f);

    % We are dealing with a real signal, so we only care about positive, real
    % freqency components.  To get the power, we want to "fold" both sides of
    % the spectrum together.  Since they are the same, we can just take one
    % side and multiply by 2.
    % To preserve variance, we have to normalize this result by dividing by
    % N^2, due to the Matlab implementation of fft.
    ehat_real = 2 * abs(a(1:length(freq_bins_real))).^2 ./ delta_f ./ ...
        (N^2);

    % Also, figure out the full spectrum
    ehat = abs(a(1:length(freq_bins)-1)).^2 ./ delta_f ./ (N^2);

    % The first frequency bin will hold the mean, not a component of the
    % variance.  Since we are interested in the variance spectrum, drop
    % the first bin.
    ehat = ehat(2:end);
    ehat_real = ehat_real(2:end);
    freq_bins = freq_bins(2:end);
    freq_bins_real = freq_bins_real(2:end);

end

function [ehat_real, ehat, freq_bins_real, freq_bins, delta_f] = ...
    varspec_est_avg(data, sample_interval, set_size, number_of_sets)

    % Get the variance spectral densitity for each subset.
    for i=1:number_of_sets
        start_idx = (i-1) * set_size + 1;
        end_idx = start_idx + set_size - 1;

        [ehat_parts_real(i, :), ehat_parts(i, :), freq_bins_real, freq_bins, delta_f] = ...
            varspec_est_real(data(start_idx:end_idx), sample_interval);
    end

    % form the average
    % Get the average values (average of columns)
    ehat = mean(ehat_parts);
    ehat_real = mean(ehat_parts_real);

end

function [ehat_real, ehat, freq_bins_real, freq_bins, delta_f] = ...
    varspec_est_multidof(data, sample_interval, dof)

    % We get 2 DOF from each set of data (thanks to real/imag components.
    num_sets = ceil(dof / 2);

    len_complete = size(data, 1);
    points_per = floor(len_complete / num_sets);

    [ehat_real, ehat, freq_bins_real, freq_bins, delta_f] = ...
        varspec_est_avg(data, sample_interval, points_per, num_sets);
end
quad =

  Columns 1 through 16

   -0.0017    0.0011   -0.0039    0.0021    0.0584    0.0736    0.0226    0.0163    0.0142    0.0103    0.0122    0.0425    0.0374    0.0372    0.0061    0.0019

  Columns 17 through 32

   -0.0021   -0.0049   -0.0024   -0.0050   -0.0068   -0.0060   -0.0011    0.0003    0.0037    0.0023    0.0016    0.0012    0.0022    0.0012    0.0002    0.0018

  Columns 33 through 48

   -0.0003   -0.0006    0.0007    0.0003   -0.0009    0.0007   -0.0012   -0.0003   -0.0002    0.0006    0.0006    0.0004    0.0002    0.0001   -0.0010   -0.0010

  Columns 49 through 60

   -0.0010    0.0023   -0.0012   -0.0020    0.0006    0.0008   -0.0009    0.0008    0.0007   -0.0004   -0.0009         0